<!DOCTYPE html>
<html>
<head><meta charset="utf-8" />
<title>solverintroduction</title><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>

<style type="text/css">
/* Overrides of notebook CSS for static HTML export */
body {
  overflow: visible;
  padding: 8px;
}
div#notebook {
  overflow: visible;
  border-top: none;
}@media print {
  div.cell {
    display: block;
    page-break-inside: avoid;
  } 
  div.output_wrapper { 
    display: block;
    page-break-inside: avoid; 
  }
  div.output { 
    display: block;
    page-break-inside: avoid; 
  }
}
</style>

<!-- Custom stylesheet, it must be in the parent directory as the html file -->
<link rel="stylesheet" type="text/css" media="all" href="../doc.css" />
<link rel="stylesheet" type="text/css" media="all" href="doc.css" />

<!-- Loading mathjax macro -->
<!-- Load mathjax -->
    <script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML"></script>
    <!-- MathJax configuration -->
    <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
        tex2jax: {
            inlineMath: [ ['$','$'], ["\\(","\\)"] ],
            displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
            processEscapes: true,
            processEnvironments: true
        },
        // Center justify equations in code and markdown cells. Elsewhere
        // we use CSS to left justify single line equations in code cells.
        displayAlign: 'center',
        "HTML-CSS": {
            styles: {'.MathJax_Display': {"margin": 0}},
            linebreaks: { automatic: true }
        }
    });
    </script>
    <!-- End of mathjax configuration --></head>
<body>
  <div tabindex="-1" id="notebook" class="border-box-sizing">
    <div class="container" id="notebook-container">

<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="Solvers-for-Linear-Algebraic-Equations">Solvers for Linear Algebraic Equations<a class="anchor-link" href="#Solvers-for-Linear-Algebraic-Equations">&#182;</a></h1>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>We have implemented multigrid solvers for linear algebraic systems arising from various finite element methods. Here we briefly present the usage for a symmetric and positive definite matrix equation <code>Ax = b</code>.</p>

<pre><code>x = A\b;
x = amg(A,b);
x = mg(A,b,elem);

</code></pre>
<p>Backslash <code>\</code> is the build-in direct solver of MATLAB. It is suitable for small size problems. <code>x = amg(A,b)</code> implements algebraic multigrid solver. To acheive multigrid efficiency, a hierarchical 'grids' is generated from the graph of <code>A</code>. When the mesh is avaiable, <code>x = mg(A,b,elem)</code> implements geometric multigrid solvers. Inside <code>mg</code>, an coarsening algorithm is applied to the mesh given by <code>elem</code>.</p>
<p>More options is allowed in <code>mg</code> and <code>amg</code>. Try <code>help mg</code> and <code>help amg</code> for possible options including tolerance, V or W-cycles, number of smoothings steps, and print level etc.</p>
<p>Here we include several examples for discrete Poisson matrices. Solvers for other finite element methods and other equations can be found</p>
<ul>
<li><a href="solverlist.html">List of Examples</a></li>
</ul>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Example:-2-D-Linear-Element">Example: 2-D Linear Element<a class="anchor-link" href="#Example:-2-D-Linear-Element">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[1]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">% setting</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">shape</span> <span class="p">=</span> <span class="s">&#39;square&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">type</span> <span class="p">=</span> <span class="s">&#39;uniform&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="nb">size</span> <span class="p">=</span> <span class="mf">2e5</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="s">&#39;Poisson&#39;</span><span class="p">;</span>
<span class="n">fem</span> <span class="p">=</span> <span class="s">&#39;P1&#39;</span><span class="p">;</span>
<span class="c">% get the matrix</span>
<span class="p">[</span><span class="n">eqn</span><span class="p">,</span><span class="n">T</span><span class="p">]</span> <span class="p">=</span> <span class="n">getfemmatrix</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">fem</span><span class="p">);</span>
<span class="c">% compare solvers</span>
<span class="n">tic</span><span class="p">;</span> <span class="nb">disp</span><span class="p">(</span><span class="s">&#39;Direct solver&#39;</span><span class="p">);</span> <span class="n">x1</span> <span class="p">=</span> <span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="o">\</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">;</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x2</span> <span class="p">=</span> <span class="n">mg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">,</span><span class="n">T</span><span class="p">.</span><span class="n">elem</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x3</span> <span class="p">=</span> <span class="n">amg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">format</span> <span class="n">shorte</span>
<span class="n">fprintf</span><span class="p">(</span><span class="s">&#39;Difference between direct and mg, amg solvers %0.2g, %0.2g \n&#39;</span><span class="p">,</span><span class="c">...</span>
         <span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x2</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">),</span><span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x3</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">));</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Direct solver
Elapsed time is 1.293569 seconds.

 Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:   263169,  #nnz:  1303561, smoothing: (1,1), iter: 10,   err = 1.35e-09,   time =  1.1 s
Elapsed time is 1.186739 seconds.

 Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method
  nnz/N: 4.96,   level:  6,   coarse grid 169,   nnz/Nc 9.38
#dof:  263169,    iter: 14,   err = 5.5672e-09,   time = 4.65 s
 
Elapsed time is 4.334722 seconds.
Difference between direct and mg, amg solvers 1.4e-09, 7.8e-08 
</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_png output_subarea ">
<img src="
B3RJTUUH4gsWBCE4Nm4s2gAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ
bmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMDozMzo1Ni9aGsoAAA1Z
SURBVHic7do9jh3HGYXhvoZyewEkwFyJUjKZCRUqMLQOAbMBbUCwl+DUjBQqnEnElAlzAjQUE9zA
dTAEf4Z37m9311d1nmcFJ6oX3VWb7XY7AUBrf2s9AACmSZAAKEKQAChBkAAoQZAAKEGQAChBkAAo
QZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChB
kAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKOG71gPGd3t7e3d313oF0MzV1dX19XXrFR3Y
bLfb1hsGt9lsrm+et15BG+/ffXj98s2z50+fvXjSegtt3P72apomJ+0xBGlxm83m179uWq+ggffv
Pvznn//9x5O/P3vx5PrmRes5NPD7L39M0/T65Rsn7THcIcEi7mv0079+9G0U675GP/37x9ZDuiFI
ML8vavS09RbaUKMzCBLMTI1Qo/MIEsxJjVCjswkSzEaNUKNLCBLMQ41QowsJEsxAjVCjywkSXEqN
UKNZCBJcRI1Qo7kIEpxPjVCjGQkSnEmNUKN5CRKcQ41Qo9kJEpxMjVCjJQgSnEaNUKOFCBKcQI1Q
o+UIEhxLjVCjRQkSHEWNUKOlCRIcpkao0QoECQ5QI9RoHYIE+6gRarQaQYJHqRFqtCZBgt3UCDVa
mSDBDmqEGq1PkOAhNUKNmhAk+IoaoUatCBJ8pkaoUUOCBB+pEWrUliDBNKkRalSAIIEaoUYlCBLp
1Ag1KkKQiKZGqFEdgkQuNUKNShEkQqkRalSNIJFIjVCjggSJOGqEGtUkSGRRI9SoLEEiiBqhRpUJ
EinUCDUqTpCIoEaoUX2CxPjUCDXqgiAxODVCjXohSIxMjVCjjggSw1Ij1KgvgsSY1Ag16o4gMSA1
Qo16JEiMRo1Qo04JEkNRI9SoX4LEONQINeqaIDEINUKNeidIjECNUKMBCBLdUyPUaAyCRN/UCDUa
hiDRMTVCjUYiSPRKjVCjwQgSXVIj1Gg8gkR/1Ag1GpIg0Rk1Qo1GJUj0RI1Qo4EJEt1QI9RobIJE
H9QINRqeINEBNUKNEggS1akRahRCkChNjVCjHIJEXWqEGkURJIpSI9QojSBRkRqhRoEEiXLUCDXK
JEjUokaoUSxBohA1Qo2SCRJVqBFqFE6QKEGNUCMEifbUCDViEiSaUyPUiHuCREtqhBrxiSDRjBqh
RnxJkGhDjVAjHhAkGlAj1IhvCRJrUyPUiJ0EiVWpEWrEYwSJ9agRasQegsRK1Ag1Yj9BYg1qhBpx
kCCxODVCjTiGILEsNUKNOJIgsSA1Qo04niCxFDVCjTiJILEINUKNOJUgMT81Qo04gyAxMzVCjTiP
IDEnNUKNOJsgMRs1Qo24hCAxDzVCjbiQIDEDNUKNuJwgcSk1Qo2YhSBxETVCjZiLIHE+NUKNmJEg
cSY1Qo2YlyBxDjVCjZidIHEyNUKNWIIgcRo1Qo1YiCBxAjVCjViOIHEsNUKNWJQgcRQ1Qo1YmiBx
mBqhRqxAkDhAjVAj1iFI7KNGqBGrESQepUaoEWsSJHZTI9SIlQkSO6gRasT6BImH1Ag1oglB4itq
hBrRiiDxmRqhRjQkSHykRqgRbQkS06RGqBEFCBJqhBpRgiClUyPUiCIEKZoaoUbUIUi51Ag1ohRB
CqVGqBHVCFIiNUKNKEiQ4qgRakRNgpRFjVAjyhKkIGqEGlGZIKVQI9SI4gQpghqhRtQnSONTI9SI
LgjS4NQINaIXgjQyNUKN6IggDUuNUCP6IkhjUiPUiO4I0oDUCDWiR4I0GjVCjeiUIA1FjVAj+iVI
41Aj1IiuCdIg1Ag1oneCNAI1Qo0YgCB1T41QI8YgSH1TI9SIYQhSx9QINWIkgtQrNUKNGIwgdUmN
UCPGI0j9USPUiCEJUmfUCDViVILUEzVCjRiYIHVDjVAjxiZIfVAj1IjhCVIH1Ag1IoEgVadGqBEh
BKk0NUKNyCFIdakRakQUQSpKjVAj0ghSRWqEGhFIkMpRI9SITIJUixqhRsQSpELUCDUimSBVoUao
EeEEqQQ1Qo1AkNpTI9QIJkFqTo1QI7gnSC2pEWoEnwhSM2qEGsGXBKkNNUKN4AFBakCNUCP4liCt
TY1QI9hJkFalRqgRPEaQ1qNGqBHsIUgrUSPUCPYTpDWoEWoEBwnSGtQo3OuXbyY1gkM22+229YbB
bTab1hOAxpy0x/iu9YAIv/5103oCbfz+yx9vX7374efvr29etN5CA/e/69+/+9B6SB/8soOl3N8b
/fDz962H0Many+PWQ7ohSLAIrxjCecp0BkGC+alRODU6jyDBzNQonBqdTZBgTmoUTo0uIUgwGzUK
p0YXEiSYhxqFU6PLCRLMQI3CqdEsBAkupUbh1GguggQXUaNwajQjQYLzqVE4NZqXIMGZ1CicGs1O
kOAcahROjZYgSHAyNQqnRgsRJDiNGoVTo+UIEpxAjcKp0aIECY6lRuHUaGmCBEdRo3BqtAJBgsPU
KJwarUOQ4AA1CqdGqxEk2EeNwqnRmgQJHqVG4dRoZYIEu6lRODVanyDBDmoUTo2aECR4SI3CqVEr
ggRfUaNwatSQIMFnahROjdoSJPhIjcKpUXOCBNOkRvHUqAJBAjVKp0ZFCBLp1CicGtUhSERTo3Bq
VIogkUuNwqlRNYJEKDUKp0YFCRKJ1CicGtUkSMRRo3BqVJYgkUWNwqlRZYJEEDUKp0bFCRIp1Cic
GtUnSERQo3Bq1AVBYnxqFE6NeiFIDE6NwqlRRwSJkalRODXqiyAxLDUKp0bdESTGpEbh1KhHgsSA
1CicGnVKkBiNGoVTo34JEkNRo3Bq1DVBYhxqFE6NeidIDEKNwqnRAASJEahRODUagyDRPTUKp0bD
ECT6pkbh1GgkgkTH1CicGg1GkOiVGoVTo/EIEl1So3BqNCRBoj9qFE6NRiVIdEaNwqnRwASJnqhR
ODUamyDRDTUKp0bDEyT6oEbh1CiBINEBNQqnRiEEierUKJwa5RAkSlOjcGoURZCoS43CqVEaQaIo
NQqnRoEEiYrUKJwaZRIkylGjcGoUS5CoRY3CqVEyQaIQNQqnRuEEiSrUKJwaIUiUoEbh1IhJkKhA
jcKpEfcEicbUKJwa8Ykg0ZIahVMjviRINKNG4dSIBwSJNtQonBrxLUGiATUKp0bsJEisTY3CqRGP
ESRWpUbh1Ig9BIn1qFE4NWI/QWIlahROjThIkFiDGoVTI44hSCxOjcKpEUcSJJalRuHUiOMJEgtS
o3BqxEkEiaWoUTg14lSCxCLUKJwacQZBYn5qFE6NOI8gMTM1CqdGnE2QmJMahVMjLiFIzEaNwqkR
FxIk5qFG4dSIywkSM1CjcGrELASJS6lRODViLoLERdQonBoxI0HifGoUTo2YlyBxJjUKp0bMTpA4
hxqFUyOWIEicTI3CqRELESROo0bh1IjlCBInUKNwasSiBIljqVE4NWJpgsRR1CicGrECQeIwNQqn
RqxDkDhAjcKpEasRJPZRo3BqxJoEiUepUTg1YmWCxG5qFE6NWJ8gsYMahVMjmhAkHlKjcGpEK4LE
V9QonBrRkCDxmRqFUyPaEiQ+UqNwakRzgsQ0qVE8NaICQUKN0qkRRQhSOjUKp0bUIUjR1CicGlGK
IOVSo3BqRDWCFEqNwqkRBQlSIjUKp0bUJEhx1CicGlGWIGVRo3BqRGWCFESNwqkRxQlSCjUKp0bU
J0gR1CicGtEFQRqfGoVTI3ohSINTo3BqREcEaWRqFE6N6IsgDUuNwqkR3RGkMalRODWiR4I0IDUK
p0Z0SpBGo0bh1Ih+CdJQ1CicGtE1QRqHGoVTI3onSINQo3BqxAAEaQRqFE6NGIMgdU+NwqkRwxCk
vqlRODViJILUMTUKp0YMRpB6pUbh1IjxCFKX1CicGjEkQeqPGoVTI0YlSJ1Ro3BqxMAEqSdqFE6N
GJsgdUONwqkRwxOkPqhRODUigSB1QI3CqREhBKk6NQqnRuQQpNLUKJwaEUWQ6lKjcGpEGkEqSo3C
qRGBBKkiNQqnRmQSpHLUKJwaEUuQalGjcGpEMkEqRI3CqRHhBKkKNQqnRiBIJahRODWCSZAqUKNw
agT3BKkxNQqnRvCJILWkRuHUCL4kSM2oUTg1ggcEqQ01CqdG8C1BakCNwqkR7CRIa1OjcGoEjxGk
ValRODWCPQRpPWoUTo1gP0FaiRqFUyM4SJDWoEbh1AiOsdlut603DG6z2UzTdH3zvPUQ2nj75//e
vnr37PnTZy+etN5CG7e/vXLSHkOQFnd7e3t3d9d6BdDM1dXV9fV16xUdECQASnCHBEAJggRACYIE
QAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRA
CYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJ
ggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACf8H1Hg8
qOD7E6gAAAAASUVORK5CYII=
"
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>For problem size of $2.6 \times 10^5$, <code>mg</code> ties with direct solver <code>\</code>. But <code>amg</code> is aroud 3-4 times slover.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Example:-2-D-Adaptive-Meshes">Example: 2-D Adaptive Meshes<a class="anchor-link" href="#Example:-2-D-Adaptive-Meshes">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[7]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">%% Lshape adaptive grids</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">shape</span> <span class="p">=</span> <span class="s">&#39;Lshape&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">type</span> <span class="p">=</span> <span class="s">&#39;adaptive&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="nb">size</span> <span class="p">=</span> <span class="mf">2e4</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="s">&#39;Poisson&#39;</span><span class="p">;</span>
<span class="n">fem</span> <span class="p">=</span> <span class="s">&#39;P1&#39;</span><span class="p">;</span>
<span class="c">% get the matrix</span>
<span class="p">[</span><span class="n">eqn</span><span class="p">,</span><span class="n">T</span><span class="p">]</span> <span class="p">=</span> <span class="n">getfemmatrix</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">fem</span><span class="p">);</span>
<span class="c">% compare solvers</span>
<span class="n">tic</span><span class="p">;</span> <span class="nb">disp</span><span class="p">(</span><span class="s">&#39;Direct solver&#39;</span><span class="p">);</span> <span class="n">x1</span> <span class="p">=</span> <span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="o">\</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">;</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x2</span> <span class="p">=</span> <span class="n">mg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">,</span><span class="n">T</span><span class="p">.</span><span class="n">elem</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x3</span> <span class="p">=</span> <span class="n">amg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">fprintf</span><span class="p">(</span><span class="s">&#39;Difference between direct and mg, amg solvers %0.2g, %0.2g \n&#39;</span><span class="p">,</span><span class="c">...</span>
         <span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x2</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">),</span><span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x3</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">));</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Direct solver
Elapsed time is 3.347438 seconds.

 Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:   738561,  #nnz:  3682043, smoothing: (1,1), iter: 10,   err = 6.29e-09,   time =  3.1 s
Elapsed time is 2.843210 seconds.

 Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method
  nnz/N: 4.99,   level:  6,   coarse grid 453,   nnz/Nc 9.72
#dof:  738561,    iter: 15,   err = 2.9982e-09,   time = 14.9 s
 
Elapsed time is 13.711571 seconds.
Difference between direct and mg, amg solvers 1.2e-08, 6e-08 
</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_png output_subarea ">
<img src="
B3RJTUUH4gsWBC4GcJctvgAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ
bmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMDo0NjowNtfTwq0AACAA
SURBVHic7d2/bmNHloDxq0XHFOOFLsBcGygzqGSUTYcN2KtOG1h0Lq8ewH4BwcqddNqEB1CobOiE
RGcMhjmBK2xM6QW0wbFpDUWR99atP6dOfb9oEjXLQg+/Jm9VnaPn5+cKAIDU/iP1AgAAqCqCBABQ
giABAFQgSAAAFQgSAEAFggQAUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQ
gSABAFQgSAAAFQgSAEAFggQAUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQ
gSABAFQgSAAAFQgSAECFd6kXYN90Ov39999TrwJAMn/7298uLi5SryIDR8/Pz6nXYNzR0dGwHpxd
nqZeSBrr5mkxWRb7G5D//Kqqzi5Ph/Ug9XISWM0eVvOm2P/8qqqmN/OqqninbYNPSDHIe/HF9Xnq
hSSwmjWreTMa11V5v4F18/Tlh6+jcT2sB+vm6eL6vLQ35dWsWc0eRuP67PJ0dF6nXk4Cd1f3Z5en
8o8SHMQzpBjkjXh6M0u9kDSGJ8cX1+fr5qmo34DU6MMv70fnJ8N68OH2/d3V/bp5Sr2ueFazZnoz
//SPy9QLSebu6r6qqg+371MvJBsEKZLSm1QPimrSixr98bGgtCZRI2rkgCDFQ5MKadLrGolymkSN
qJEbghQVTTLfpLdqJEpoEjWiRs4IUmw0yXCT9tdI2G4SNaJGfRCkBGiSySa1qZGw2iRqRI16Ikhp
0CRjTWpfI2GvSdSIGvVHkJKhSWaa1LVGwlKTqBE18oIgpUSTDDTJrUbCRpOoETXyhSAlRpOyblKf
Goncm0SNqJFHBCk9mpRpk/rXSOTbJGpEjfwiSCrQpOya5KtGIscmUSNq5B1B0oImZdQkvzUSeTWJ
GlGjEAiSIjQpiyaFqJHIpUnUiBoFQpB0oUnKmxSuRkJ/k6gRNQqHIKlDk9Q2KXSNhOYmUSNqFBRB
0ogmKWxSnBoJnU2iRtQoNIKkFE1S1aSYNRLamkSNqFEEBEkvmqSkSfFrJPQ0iRpRozgIkmo0KXmT
UtVIaGgSNaJG0RAk7WhSwialrZFI2yRqRI1iIkgZoElJmqShRiJVk6gRNYqMIOWBJkVukp4aifhN
okbUKD6ClA2aFK1J2mokYjaJGlGjJAhSTmhShCbprJGI0yRqRI1SIUiZoUlBm6S5RiJ0k6gRNUqI
IOWHJgVqkv4aiXBNokbUKC2ClCWa5L1JudRIhGgSNaJGyRGkXNEkj03Kq0bCb5OoETXSgCBljCZ5
aVKONRK+mkSNqJESBClvNKlnk/KtkejfJGpEjfQgSNmjSc5Nyr1Gok+TqBE1UoUgWUCTHJpko0bC
rUnUiBppQ5CMoEmdmmSpRqJrk6gRNVKIINlBk1o2yV6NRPsmUSNqpBNBMoUmHWyS1RqJNk2iRtRI
LYJkDU3a0yTbNRL7m0SNqJFmBMkgmrSzSSXUSLzVJGpEjZQjSDbRpK0mlVMj8bpJ1Iga6UeQzKJJ
myaVViPxsknUiBplgSBZRpOkSbff/VpajcSmSdSookY5IEjGFd6kqqoWk+WwHqzmTeqFpLFuHlfz
Zv3wGG32uTbUKCMEyb5imyTf1J1dng5PjuPMPtdGvqk7uzw9uzyNM/tcG2qUF4JUhAKb9PK5UZzZ
59q8fG40rI8jzD7XhhplhyCVoqgmvd7FUFqTXu9iCD37XBtqlCOCVJBCmvTWnrpymvTWnrpymkSN
MkWQymK+Sft3eJfQpP07vEtoEjXKF0EqjuEmtTlvZLtJbc4b2W4SNcoaQSqRySa1P/1qtUntT79a
bRI1yh1BKpSxJnW9i8Fek7rexWCvSdTIAIJULjNNcrsZyFKT3G4GstQkamQDQSqagSb1uafORpP6
3FNno0nUyAyCVLqsm9T/1tTcm9T/1tTcm0SNLCFIyLVJvu7wzrdJvu7wzrdJ1MgYgoSqyrBJfidK
5NgkvxMlcmwSNbKHIOEPGTUpxHyjvJoUYr5RXk2iRiYRJPwliyaFm7aXS5PCTdvLpUnUyCqChH+j
vEmhZ7/qb1Lo2a/6m0SNDCNI2Ka2SXEmkWtuUpxJ5JqbRI1sI0jYQWGT4tRI6GxSnBoJnU2iRuYR
JOymqkkxayS0NSlmjYS2JlGjEhAkvElJk+LXSOhpUvwaCT1NokaFIEjYJ3mTUtVIaGhSqhoJDU2i
RuUgSDggYZPS1kikbVLaGom0TaJGRSFIOCxJkzTUSKRqkoYaiVRNokalIUhoJXKT9NRIxG+SnhqJ
+E2iRgUiSGgrWpO01UjEbJK2GomYTaJGZSJI6CBCk3TWSMRpks4aiThNokbFIkjoJmiTNNdIhG6S
5hqJ0E2iRiUjSOgsUJP010iEa5L+GolwTaJGhSNIcOG9SbnUSIRoUi41EiGaRI1AkODIY5PyqpHw
26S8aiT8NokaoSJI6MNLk3KskfDVpBxrJHw1iRpBECT00rNJ+dZI9G9SvjUS/ZtEjbBBkNCXc5Ny
r5Ho06TcayT6NIka4SWCBA8cmmSjRsKtSTZqJNyaRI2whSDBj05NslQj0bVJlmokujaJGuE1ggRv
WjbJXo1E+ybZq5Fo3yRqhJ0IEnw62CSrNRJtmmS1RqJNk6gR3kKQ4NmeJtmukdjfJNs1EvubRI2w
B0GCfzubVEKNxFtNKqFG4q0mUSPsR5AQxFaTyqmReN2kcmokXjeJGuEggoRQ/mzSvKqqomokNk1a
TJar2UNRNRJ/NenhUf4aUCPsd/T8/Jx6DcYdHR2NxgW9EW9ZzZuqqob1YHhynGQB64fHqqoSvrp8
REj4d4DfQHKrecM7bRvvUi+gCKt58+m3sv51LNYPT/J+NDw5vrgeJ1nDat4sJsuqqpIsYDVv5MNB
qgVUVTW9mY/OT1L1YHozlyCdfTwdngySrCGtux/vUy8hGwQphg+37+9+vL/69jn1QqJaN093P96f
XZ6uZg+j85PVvJEv8WIv4+FpeHKcZAGrWbOaPVxcj+UdOdVvYDhZDuvjJN+X3l3dD+vBqKrPPp4u
vi4/3L4f1mU16cv3kw+/vP/ywyT1QvLAM6QYzi5PL67Pb7/7NfVC4vlrF8O4rqLMPt9jWA/OLv8r
8gJe7mKIM/tcm5e7GIYnMWafa/Pl+8nF9bioR6c9EaRIimrSzj11RTXp9Z660pr0ek9d6Nnn2lAj
BwQpnkKatGeHdyFNemuHdzlNemuHdzlNokZuCFJU5pt08LyR+SbtP29UQpP2nzcqoUnUyBlBis1w
k1qefjXcpDanX203qc3pV9tNokZ9EKQETDap010MJpvU/i4Gq01qfxeD1SZRo54IUhrGmuRwM5Cx
JnW9Gchek7reDGSvSdSoP4KUjJkmOd9TZ6ZJbvfUWWqS2z11lppEjbwgSCkZaFLPW1MNNKnPrak2
mtTn1lQbTaJGvhCkxLJukpc7vLNuUv87vHNvUv87vHNvEjXyiCCll2mTPE6UyLRJviZK5NskXxMl
8m0SNfKLIKmQXZO8zzfKrkl+5xvl2CS/841ybBI18o4gaZFRkwJN28uoSSGm7eXVpBDT9vJqEjUK
gSApkkWTgs5+zaJJ4Wa/5tKkcLNfc2kSNQqEIOmivEkRJpErb1LoSeT6mxR6Ern+JlGjcAiSOmqb
FKFGQm2TQtdoswC1TQpdI6G5SdQoKIKkkcImRauRUNikODXaLEBhk+LUSOhsEjUKjSAppapJkWsk
VDUpZo02C1DVpJg1EtqaRI0iIEh6KWlSkhoJJU26u7qPXKPNApQ0KX6NhJ4mUaM4CJJqyZuUsEYi
eZOG9fFishydn6RaQPImpaqR0NAkahQNQdIuYZOS10gkbNJq1iy+Lof1INUCqtRNSlsjkbZJ1Cgm
gpSBJE1SUiORpEkvnxvFmX3+llRN0lAjkapJ1CgygpSHyE1SVSMRuUlbuxgizD7fL36T9NRIxG8S
NYqPIGUjWpMU1khEa9LOPXVFNUlbjUTMJlGjJAhSTiI0SW2NRIQm7dnhXUiTdNZIxGkSNUqFIGUm
aJOU10gEbdLB80bmm6S5RiJ0k6hRQgQpP4GalEWNRKAmtTz9arhJ+mskwjWJGqVFkLLkvUkZ1Uh4
b1KnuxhMNimXGokQTaJGyRGkXHlsUnY1Eh6b5HAzkLEm5VUj4bdJ1EgDgpQxL03KtEbCS5Oc76kz
06QcayR8NYkaKUGQ8tazSVnXSPRsUs9bUw00Kd8aif5NokZ6EKTsOTfJQI2Ec5O83OGddZNyr5Ho
0yRqpApBssChSWZqJBya5HGiRKZNslEj4dYkaqQNQTKiU5OM1Uh0apL3+UbZNclSjUTXJlEjhQiS
HS2bZLJGomWTAk3by6hJ9mok2jeJGulEkEw52CTDNRIHmxR09msWTbJaI9GmSdRILYJkzZ4mma+R
2NOkCJPIlTfJdo3E/iZRI80IkkE7m1RIjcTOJkWokVDbpBJqJN5qEjVSjiDZtNWkomoktpoUrUZC
YZPKqZF43SRqpB9BMmvTpAJrJKRJi6/L1TxqjcSmSYvJMubrvlyANGk1byRL5dRIvGwSNcrC0fPz
c+o1GHd0dPTz/12nevXpzWx6Mx/Wg1QLWDdPCV9dFlBVVZ819PlP6P/qPW0+HxT7d0B+Ax9u359d
nqZaw8//ecM7bRvvUi8AAa2bp8VkeXZ5upgsr759TrSGxy8/TEbjOsk/z1fz5u7qflgPRuNaPjA5
+PLD13Xz9Om3y2F93PVnF5N/TW/mVVX1WUAf05uZfERLtYCqqu6u7s8+no7GaT6dyKsvvi5H4zrt
v41wEEEy6+U3daPz+ssPX5M0ad08jsb16PxkejOL3KTVrFl8XX64fb+aNcN6sJj8y/kd+erb5y8/
fP3028dO72iLyXI1e7i4HstHhD4LcCPPjc4uT4f1YN08xV/AxvBkkCQGm2/qRuP67ur+w+17mqQZ
z5Bs2npuFGH2+X4X1+fDeiDvj3Fs7WLoeQfrsB58+u2jfFRq+SOLyXLxdSkLSLLH4eUuhmF9HGH2
uTYvnxvFmX2OngiSQTt3MRTVpJ176mI26WWNNj8es0mv99SFnn2uzetdDDRJP4JkzZ49dYU0ac8O
7zhNel2jzY/HadJbO7zLadJbe+poknIEyZSDO7zNN+ngeaPQTXqrRpsfD92k/eeNSmjS/h3eNEkz
gmRHy/NGhpvU8vRruCbtr9Hmx8M1qc3pV9tNanPeiCapRZCM6HT61WSTOt3FEKJJbWq0+fEQTWp/
F4PVJrU//UqTdCJIFjjcxWCsSQ43A/ltUvsabX7cb5O63gxkr0ld72KgSQoRpOw53wxkpknO99T5
atL0ZtapRpsf99Ukt3vqLDXJ7WYgmqQNQcpbz3vqDDSp562p/Zt0dnk6vZm7nfn10qQ+t6baaFKf
e+pokioEKWNebk3Nukle7vDu0yS5i0HucXB7R+vZpP53eOfepP63ptIkPQhSrjze4Z1pkzxOlHBr
0ua5kcM9Di85N8nXRIl8m+TrDm+apARBypL3iRLZNcn7fKOuTdraxRC/SX7nG+XYJL8TJWiSBgQp
P4HmG2XUpEDT9to3aeeeuphNCjFtL68mhZhvRJOSI0iZCTptL4smBZ392qZJe3Z4x2lSuNmvuTQp
3LQ9mpQWQcpJhNmvypsUYRL5/iYdPG8UukmhJ5Hrb1Lo2a80KSGClI1ok8jVNilCjTYLqHYloeXp
13BNCl2jzQLUNinOJHKalApBykO0GgmFTYpWo80Cqn9PQqe7GEI0KU6NNgtQ2KQ4NRI0KQmClIHI
NRKqmhS5RpsFVH8moevNQJXvJsWs0WYBqpoUs0aCJsXHCHPtktRInF2eVlV1+92vSWafV1V1cX0+
vZl9+X6yfni8uD5fTJYOf8hq1qybJ7efHdbH05uZvCm7LUDucTi7PHVbQPXn2Vt5c3T7E5xJk6Y3
M/nPj/zqL8Wvkdg0idnncRAk1RLWSCRv0mhcy1v5unl0+xPWzdNq3gzrQc83FLcFbK5erarKbQHr
5mndPF2cjx1+tj8NTUpVI0GTYiJIeiWvkUjYJPmm7urbZ/mY4vYRYVFLz57OLk+7/iYXk+Xw5Lj6
85fQ9R153TzdfvfraFwP68Fq3ny6/tj1He3u6n7z6SpVEtI2KW2NBE2KhmdISimpkUjyPOnlc6Oe
d7DKG8r0Zr6aNe1/6uVzI4e7haRGn367HJ2fuD1P2jw3ijb7/C2pnidpqJHgeVIcBEkjVTUSkZv0
ehdD5Ca93sXQqUkvalRvFtCpSVu7GApskp4aCZoUAUFSR2GNRLQmvbWnLlqT3tpT17JJr2u0WUDL
Ju3cU1dUk7TVSNCk0AiSLmprJCI0af8O7whN2r/D+2CT3qrRZgEHm7Rnh3chTdJZI0GTgiJIiiiv
kQjapDbnjYI2qc15oz1N2l+jzQL2NOngeSPzTdJcI0GTwiFIWmRRIxGoSe1PvwZqUvvTrzub1KZG
mwXsbFLL06+Gm6S/RoImBUKQVMioRsJ7k7rexeC9SV3vYthqUvsabRaw1aROdzGYbFIuNRI0KQSC
lF52NRIem+R2M5DHJjncDFS9aFLXGm0WsGmSw81AxpqUV40ETfKOICWWaY2Elyb1uafOV5McarRZ
gFuNNgv49NtH+QU6HPs106QcayRokl8EKaWsayR6Nqn/rak9m7SaN5uVOPy43JJ3dnm6+XO6mt7M
5MejzT73q3+T8q2RoEkeEaRkDNRIODfJ1x3ezk3afFPncI9D9eK5kXy4cXhH3nxTF232eQh9mpR7
jQRN8oUgpWGmRsKhSX4nSjg06eVzI4e7hba+qXO4W+jlc6M4s8/DcWuSjRoJmuQFQUrAWI1EpyaF
mG/UqUmvdzF0atLO50admvR6F0NpTbJUI0GT+iNIsZmskWjZpHDT9lo26a09dS2btGcXQ8smvbWn
rpwm2auRoEk9EaSoDNdIHGxS6NmvB5u0f4f3wSYd3FN3sEn7d3iX0CSrNRI0qQ+CFI/5Gok9TYoz
iXxPk9qcN9rTpJY7vPc0qc15I9tNsl0jQZOcEaRICqmR2NmkODUSO5vU/vTrziZ1Om+0s0ntT79a
bVIJNRI0yQ1BiqGoGomtJsWskdhqUte7GLaa5HD6datJXe9isNekcmokaJKDo+fn59RrMO7o6GhY
D4qq0cZisry7uh/Wg+HJccwabUxvZjL/220BcqmPHHp1u4ths4DRuHa4i0H+KVNVldv8+HXztJj8
azFZXlyfyxT2yCRIi8my2P8LbP4K8U7bBkEK7ujoKPUSAKR0cXHxz3/+M/UqMvAu9QKKIP86HtaD
1AuJbTVr7n68XzdPw3rg9m/8/haT5fRm5vYBRdx+9+vw5Njt6yb51/HV//zvTz/95PbqQDl4hhRD
mV8ly3OjD7+8H43rOLPP3zIa133uu6uqyu1uIXnFJN+VATkiSDEU+HhzaxdDhNnn+0WYfb7FYaIE
UDiCFElRTdq5p66oJlEjwAFBiqeQJu3Z4V1Ik6gR4IYgRWW+SQfPG5lvEjUCnBGk2Aw3qeXpV8NN
okZAHwQpAZNN6nQXg8kmUSOgJ4KUhrEmOdwMZKxJ1AjojyAlY6ZJzvfUmWnSl+8nFTUCeiNIKRlo
Us9bUw00aVgPVvOG069AfwQpsayb5OUO76ybJD919e2zwz0OALYQpPQybZLHiRKZNmnz3MjhHgcA
rxEkFbJrkvf5Rtk1aWsXA00C+iNIWmTUpEDT9jJq0s49dTQJ6IkgKZJFk4LOfs2iSXt2eNMkoA+C
pIvyJkWYRK68SQfPG9EkwBlBUkdtkyLUSKhtUsvTrzQJcEOQNFLYpGg1Egqb1OkuBpoEOCBISqlq
UuQaCVVNcrgZ6EWTHkItEbDlXeoF4E2bJslJl1TLSFIjIdcf3H7369W3z/Ffvaqqi+vz6c1sPXla
NMuL6/H0Ztb1TxjWg8VkWf09xOoAawiSasmblLBGInmT1s3TaFyvHx6H9bHDj7v9FFAmgqRdwiYl
r5FI2CT5pu7TPy6nN7PVrHG7PnXdPPpeF2ATz5AykOR5kpIaiSTPk14+N+p5ByuANghSHiI3SVWN
ROQmvd7FQJOA0AhSNqI1SWGNRLQmvbWnjiYBQRGknERoktoaiQhN2r/DmyYB4RCkzARtkvIaiaBN
anPeiCYBgRCk/ARqUhY1EoGa1P70K00CQiBIWfLepIxqJLw3qetdDDQJ8I4g5cpjk7KrkfDYJIeb
gSqaBPhGkDLmpUmZ1kh4aZJbjQRNAjwiSHnr2aSsayR6NqlPjQRNAnwhSNlzbpKBGgnnJvWvkaBJ
gBcEyQKHJpmpkXBokq8aCZoE9EeQjOjUJGM1Ep2a5LdGgiYBPREkO1o2yWSNRMsmhaiRoElAHwTJ
lINNMlwjcbBJ4WokaBLgjCBZs6dJ5msk9jQpdI0ETQLcECSDdjapkBqJnU2KUyNBkwAHBMmmrSYV
VSOx1aSYNRI0CeiKEeZmbZp09vF0ejP79NvHmANnN9YPT1VVJXnp0bg+uzyVJIzG9Yfb95GXcXb5
X9Ob2fRmfvHz32O+LpApgmTZsB6cfTy9u7of1oMvP3xNtYx185Tw1f9Yw8Nj8jUA2I8gWbaaNYuv
y6tvn++u7j/cvh/WgyRrSPht4d3V/apq5H9fffucZA3Tm1mS1wWywzMkszYliDb7XBv5su7i+nw0
ruPMPgfQB0GyaetzSYFN2trFEGH2OYCeCJJBO78lK6pJO/fU0SRAOYJkzZ5nNoU0ac8Ob5oEaEaQ
TDm4g8B8kw6eN6JJgFoEyY6W+9kMN6nl6VeaBOhEkIzotLvaZJM63cVAkwCFCJIFDmd9jDXJ4WYg
mgRoQ5Cy53zy1EyTnO+po0mAKgQpbz3vQTDQpJ63ptIkQA+ClDEvt/Jk3SQvd3jTJEAJgpQrj3fE
ZdokjxMlaBKgAUHKkvcbS7Nrkvf5RjQJSI4g5SfQ/dkZNSnQtD2aBKRFkDITdJpDFk0KOvuVJgEJ
EaScRJgtpLxJESaR0yQgFYKUjWiT7tQ2KUKNBE0CkiBIeYg8d1Vhk6LVSNAkID6ClIEkU8BVNSly
jQRNAiIjSNolqZFQ0qQkNRI0CYiJIKmWsEYieZMS1kjQJCAagqRX8hqJhE1KXiNBk4A4CJJSSmok
kjRJSY0ETQIiIEgaqaqRiNwkVTUSNAkIjSCpo7BGIlqTFNZI0CQgKIKki9oaiQhNUlsjQZOAcAiS
IsprJII2SXmNBE0CAiFIWmRRIxGoSVnUSNAkIASCpEJGNRLem5RRjQRNArwjSOllVyPhsUnZ1UjQ
JMAvgpRYpjUSXpqUaY0ETQI8IkgpZV0j0bNJWddI0CTAF4KUjIEaCecmGaiRoEmAFwQpDTM1Eg5N
MlMjQZOA/ghSAsZqJDo1yViNBE0CeiJIsZmskWjZJJM1EjQJ6IMgRWW4RuJgkwzXSNAkwBlBisd8
jcSeJpmvkaBJgBuCFEkhNRI7m1RIjQRNAhwQpBiKqpH4q0kPT1VhNRI0CejqXeoFFOHux/sPv7xf
zZrUC4nt7OOppGhU1RfX4yS/gXXzuG6ekrz08GQwGtfTm/nFz3+P/+pAdghSDMOT4+nNPPUqUlo/
PKb6DawfHquqKvz3D2SBIMVQ1Jd1L91d3Y+qev3wODw5/nD7flgP4q9hMVmuZk2qbwu/fD8Zjesk
Lw1kh2dICEW+rLu4HkuNIsw+1+bL95OL6/Ho/CT1QoA8ECQEsbWLIcLsc23+rBEfj4C2CBL827mn
rqgmUSPAAUGCZ3t2eBfSJGoEuCFI8OngeSPzTaJGgDOCBG9ann413CRqBPRBkOBHp7sYTDaJGgE9
ESR44HAzkLEmUSOgP4KEvpzvqTPTJGoEeEGQ0EvPW1MNNIkaAb4QJLjzcod31k2iRoBHBAmOPE6U
yLRJ1AjwiyDBhff5Rtk1iRoB3hEkdBZo2l5GTaJGQAgECd0Enf2aRZOoERAIQUIHESaRK28SNQLC
IUhoK0KNhNomUSMgKIKEVqLVSChsEjUCQiNIOCxyjYSqJlEjIAKChAOS1EgoaRI1AuIgSNgnYY1E
8iZRIyAagoQ3Ja+RSNgkagTERJCwm5IaiSRNokZAZAQJO6iqkYjcJGoExEeQsE1hjUS0JlEjIAmC
hH+jtkYiQpOoEZAKQcJflNdIBG0SNQISIkj4QxY1EoGaRI2AtAgSqiqrGgnvTaJGQHIECfnVSHhs
EjUCNCBIpcu0RsJLk6gRoARBKlrWNRI9m0SNAD0IUrkM1Eg4N4kaAaoQpEKZqZFwaBI1ArQhSCUy
ViPRqUnUCFCIIBXHZI1EyyZRI0AnglQWwzUSB5tEjQC1CFJBzNdI7GkSNQI0I0ilKKRGYmeTqBGg
HEEqQlE1Ei+a9FhRIyAHBMm+AmskpEmLyXIxWVIjQL93qReAsIqtkRjWg+HJ8fDkeFgfp14LgAP4
hGRZ4TWq/vymLubscwDOCJJZ1Gjz3Cja7HMAfRAkm6jR1i4GmgToR5AMokY799TRJEA5gmQNNdqz
w5smAZoRJFOo0cHzRjQJUIsg2UGNWp5+pUmATgTJCGrU6S4GmgQoRJAsoEYONwPRJEAbgpQ9auR8
Tx1NAlQhSHmjRj1vTaVJgB4EKWPUyMsd3jQJUIIg5YoaeZwoQZMADQhSlqiR9/lGNAlIjiDlhxoF
mrZHk4C0CFJmqFHQ2a80CUiIIOWEGkWYRE6TgFQIUjaoUYQaCZoEJEGQ8kCNotVI0CQgPoKUAWoU
uUaCJgGRESTtqFGSGgmaBMREkFSjRglrJGgSEA1B0osaJa+RoElAHARJKWqkpEaCJgERECSNqJGq
GgmaBIRGkNShRgprJGgSEBRB0oUaqa2RoElAOARJEWqkvEaCJgGBECQtqFEWNRI0CQiBIKlAjTKq
kaBJgHcEKT1qlF2NBE0C/CJIiVGjTGskaBLgEUFKiRplXSNBkwBfCFIy1MhArdmuMwAAAXpJREFU
jQRNArwgSGlQIzM1EjQJ6I8gJUCNjNVI0CSgJ4IUGzUyWSNBk4A+CFJU1MhwjQRNApwRpHiokfka
CZoEuCFIkVCjQmokaBLggCDFQI2KqpGgSUBX71IvoAiLyfLiejy9maVeSALr5mk1b0bjejVvVvMm
9XJiG9aDxWT54ef/Tr0QIANHz8/Pqddg3HQ6/f3331OvAin99NNPqZcAZIAgAQBU4BkSAEAFggQA
UIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA
UIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA
UIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA
UIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAV/h90BlRfclZCBAAAAABJRU5ErkJggg==
"
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>The finest mesh is several uniform refinement of an adaptive mesh shown above. Now the multigrid outperforms the direct solver around the size of 7e5. amg is 4-5 times slower.</p>

</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="Example:-3-D-Linear-Element">Example: 3-D Linear Element<a class="anchor-link" href="#Example:-3-D-Linear-Element">&#182;</a></h2>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In&nbsp;[10]:</div>
<div class="inner_cell">
    <div class="input_area">
<div class=" highlight hl-matlab"><pre><span></span><span class="c">% Cube uniform grids</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">shape</span> <span class="p">=</span> <span class="s">&#39;cube&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="n">type</span> <span class="p">=</span> <span class="s">&#39;uniform&#39;</span><span class="p">;</span>
<span class="n">mesh</span><span class="p">.</span><span class="nb">size</span> <span class="p">=</span> <span class="mf">2e4</span><span class="p">;</span>
<span class="n">pde</span> <span class="p">=</span> <span class="s">&#39;Poisson&#39;</span><span class="p">;</span>
<span class="n">fem</span> <span class="p">=</span> <span class="s">&#39;P1&#39;</span><span class="p">;</span>
<span class="c">% get the matrix</span>
<span class="p">[</span><span class="n">eqn</span><span class="p">,</span><span class="n">T</span><span class="p">]</span> <span class="p">=</span> <span class="n">getfemmatrix3</span><span class="p">(</span><span class="n">mesh</span><span class="p">,</span><span class="n">pde</span><span class="p">,</span><span class="n">fem</span><span class="p">);</span>
<span class="c">% compare solvers</span>
<span class="n">tic</span><span class="p">;</span> <span class="nb">disp</span><span class="p">(</span><span class="s">&#39;Direct solver&#39;</span><span class="p">);</span> <span class="n">x1</span> <span class="p">=</span> <span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="o">\</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">;</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x2</span> <span class="p">=</span> <span class="n">mg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">,</span><span class="n">T</span><span class="p">.</span><span class="n">elem</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">tic</span><span class="p">;</span> <span class="n">x3</span> <span class="p">=</span> <span class="n">amg</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">A</span><span class="p">,</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">);</span> <span class="n">toc</span><span class="p">;</span>
<span class="n">fprintf</span><span class="p">(</span><span class="s">&#39;Difference between direct and mg, amg solvers %0.2g, %0.2g \n&#39;</span><span class="p">,</span><span class="c">...</span>
         <span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x2</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">),</span><span class="n">norm</span><span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">x3</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">eqn</span><span class="p">.</span><span class="n">b</span><span class="p">));</span>
</pre></div>

</div>
</div>
</div>

<div class="output_wrapper">
<div class="output">


<div class="output_area">

<div class="prompt"></div>


<div class="output_subarea output_stream output_stdout output_text">
<pre>Direct solver
Elapsed time is 0.443978 seconds.

 Multigrid V-cycle Preconditioner with Conjugate Gradient Method
#dof:    35937,  #nnz:   202771, smoothing: (1,1), iter: 11,   err = 5.50e-09,   time = 0.18 s
Elapsed time is 0.150138 seconds.

 Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method
  nnz/N: 5.81,   level:  4,   coarse grid 305,   nnz/Nc 30.80
#dof:   35937,    iter: 10,   err = 3.5316e-09,   time = 0.92 s
 
Elapsed time is 0.645183 seconds.
Difference between direct and mg, amg solvers 1e-09, 4.2e-09 
</pre>
</div>
</div>

<div class="output_area">

<div class="prompt"></div>




<div class="output_png output_subarea ">
<img src="
B3RJTUUH4gsWBQUo2/O3LwAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ
bmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMTowNTo0MBueoAYAACAA
SURBVHic7d3Nkd02t4VhtusLoD1zBCoF0WGoetxheOChBwpDY1WHoRRcpVIEKk+sDPoOcE1DJA8O
iL+9Nvg+dQe++qRunD8sLhCHfHh7e1sAALD2i/UAAABYFgIJACCCQAIASCCQAAASCCQAgAQCCQAg
gUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAAABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAg
gUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAAABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAg
4X/WAwBw7OHhYf3vt7c3w5EAY9CQAC0P//r4/Y9lWV5enz9+/yMOJ2BWNCRAwho5IYeWZfn65dv6
H2smUZUwMQIJMBaSZs2h4OuXb++f3sV/Ev7Cw8MDmYRZEUiAjX0luuX907s1n6hKmBiBBIx2WIli
IX5+/+3Pj9//WBfuVlQlzIpAAgbJrET7xbrl55K0/hCqEiZDIAF95S/NnUVVwmQIJKCXu0tze2sN
Cut18f+0L0kBVQnTIJCAxoor0WHe5Fir0kIswTMCCWimoBKdcqskBazgwTsCCajV5CxRcT3aYAUP
fhFIQKGGuxU2abQ/gbRKl6R4PFQluEMgAaf1XpprgqoEdwgkIFenDdxnF+tySlJAVYIvBBJwX79K
1OrUUQJVCV4QSMBN/b7TmpA4gbTKL0kBVQkuEEjA1rAcGlCPYlQliCOQgP+M3K1QmUZnS1JAVYIy
AgkwWJob3I02qErQRCDhukxOEaXlnEBalZWkgKoEQQQSrsj2i0S29ShGVYIUAgkXolCJ2qZRTUkK
uDArdBBIuAQX11YwxAoeFBBImJlCJYql28ypE0ir+pK0YgUPtggkTEgthwKdU0cJVCUYIpAwFZbm
mqAqwQSBhBloVqLY3XpUtl4XNFy1W1GVMB6BBN9cVCIXi3WHqEoYiUCCS/qVaLAeJSmgKmEYAgme
eMwhv/UoRlXCAAQSfHCxNLeXmUY1J5BW/UpSQFVCbwQSpHmsRHOjKqEfAgminFai2PjFut4lKaAq
oZNfrAcA/OThXy+vz8uy/P7bn7//9qf1oErMceooIVSltcIC9WhIkLBZmvv65dunD59fXp/DnB4y
yXVbSmhyAmk1piQFVCW0RSDBWM7SXPhf207cXU1fj2KcVUIrBBJsJHYr3JrNP37/w0VVMk+jkSUp
4B4WaIJAwmjpSpSeSd1VpUthBQ+V2NSAQdbdCh+//3E3TtKRE6qS5maHs9WkU7iGktT8x+ZgswOK
0ZDQ16kvEuXP5ppVyXyxTgRVCWUIJPRy9otEBbO5l7NKJsafSdpgswPOIpDQWP21FU71Hp2qRD3a
oyrhFAIJzdRcW6FyNjevSmXj752j5iVpRVVCDgIJteorUZNJU6cqYbUeJfz+258vr89UJaQRSCik
edlTk6ok0kIOmZSkdQPk5lXgrBLSCCSc1vayp5vpsr7iDK5KxdP9fGXuVg4FIRo5q4QEAgm5elSi
fgfv5meVdAwoSemnOqzXxX9CVcIhAgl3aC7N5RhQlZQX6wZIV6K9OBqpStgjkHBT7zsSjZnN+1Ul
R2nUvCS1ekqpSogRSNgaU4kO58dObWatSku3U19nOT2BdLYSLf8+0vU6RvtopCphRSDhPxPcpDXh
4vvCa0pSQQ6dRVXCQiBhsThLZLjY1WQFz9FiXaUeq523opF7WIBAui6r3Qrms3llVTIf/wADKtEt
rOBdGYF0RZpLc4MX0wz3hVstG+as2rV9Tm490rsjYQXvmgikC1HYwC1VLwqqktT4GzKsRLdQlS6I
QLoEkUqkOZvnVyXN8Z+yqSaGOZS5yYKqdCkE0swUKtFKeTa/4AY8R5exoCpdB4E0IakcyqQQBumq
1CpQbR/pyEp095Ge2olOVboCAmkqIktze8r1KHarKnkZf0KctevXVH2hKk2PQJqBeCVyN5vPdGHW
w0qkc+O+gpFQlSZGIPkmW4lOUViv24irksjcfYrgrrmGqEqzIpBcEq9EMY+z+ap5VRoQvZkD7l2S
8h/p+6d3P/758fjr49lfQVWaD4HkiaMcClynUfDy+vz+6Z1gh9twXYm+//X349PpQFqoStMhkHyY
Y2nOnfjmPbJnlYoHpnMmqRJVaRoEkjR3lSiWOdnJlo/DuyTUxFLbR+q6Em3URyNVaQ4EkijvlWiO
Q+8986/Qts2hTiXJ6vmhKnlHIEmTbQ/Tu3vpz/EreLJrhk20ikbuYeEagSQtTHzu5iDv9SjzGmvL
kCOGTx8+x7+xuWnOJMVYwXOKQFKnfDr90KnZzWPcxvJfnYJHGn5y2OZXPMLxil/T5tHICp47BJKu
tR6Zn7e4lIILByxNX53NWaIBl/mZsiQFVCVfCCQ3XFQl7/Na8fjrX52Zds0V6xSNVCUvCCRpm3NI
4lXJexpVKn51FI4zJi5JAVXJBQLJHxdVKYdasrba5bV/dQ4f6ZSVqP417RqNVCVxBJK6w412glXJ
+/F1w/GnXx3ZHJq+JAVUJWUEkmM6VekKE9lZ+1dH5MXSNyAaqUqaCCQHEt9GEqxKHnWa/uKrDbl4
jS5SkgKqkiACaQa2ValsCtOZoPtNwevq3EU0fE2HRSNVSQqB5MPdSzbUX/qzjPcD6k5Xcgv/semv
OhmccKmSFFCVdBBIU3E08c3n7m4FnXN+LgyORqqSAgLJjfzr2g2b+IrnC5HIbDXf5T/bVkX2lAuW
pIALs5ojkOY0oCp5n7Pqx5+zgfvwJZiyyPZ4OCbRyAqeIQLJk7MX/2aNqJMmzyqvjjJW8EwQSJPr
dDB+zXrU/DutslVJZ9XOcCRUpfEIJGfK7pDU9mC8coIwn3/Pjr84h9TO+aEAVWkkAukqZA/GlQ3L
CcFX51Q16Tpy87pGVRqGQPKnrCTF/3apmGRFVnKK5Yzf6nJzVCVlVKUBCKTLqTkYnz6NzPNAqiqZ
VxOdkVztohtWCCSXakpS/BOWM5NvkxlBZKrd6FGJah6pr6qk+Zq2Er8QxFJvBNJ1SR2MD7APVNk7
QSzRq/Py+mw4DPNqEo/kxz8/Hn99HPPr9u+N8Ekhk7oikLyqL0nxz1nuTcois1Kxzfi99A9fVam3
73/9/fjUPZAOn/DrHLfZIpBwvyp5T6OVciW65eX1+f3TO8MJ8SIlyeN7Yz4EkmMfv//R9j6nvQ/G
rWbVr1++ffrweRk41zR/pLJVafBr2rwkFV//CT0QSL59+vC54Ufl8NKfIkfHZWY67DW8MKtOSWoo
85kkjUYikLAVr+A5nYbiHHL6EG652laUjfponOkYZT4Ekm+ttjbc+smLt8/tNPXu65dv8X9vHsX4
V2eCklTwjF02+K0QSLip7en0rp/tw8Neqwm05pGuORRGvv6/4T/ih6NQlUx++9lo7H0pQjREILnX
qSStn3nxqqQ8tlP2kRNegvjPbauSr5LE0pxHBBLuUzgY37g73XiZOjeVKP7z/fg3+RQIvjq9paOx
SUJf6vnUQSDNoHlJOvy01xyMN1z3WweT+Gsu0mgfLbfE869tVVIuSQ0rEWlkhUDCVmLGMTwY97I0
d/fJuVWJNn8nPekrVCXbWXuNRpbmZkIgTaLfdrtbv2sZMgUUTDeah/A5ObT+zVuLdfGfW1UlnZLU
6cvO1CNDBBJ+kjnXDDgYL5tVRebKWP7S3HJ+/ApVabDNl8ya//ApnzQvCKR51Jeks7Nh5sH4qVGF
H/jy+txjuuktfqT5lSjTrWoyviqtKTiYl2VbFCOQUKXVwXiTMwEi9ehUJdr8w+LxD65Kba9ZlZZ4
b7RdP6QemSOQplJTkmo+2JUb8Mr+4YZ5GoU8KB5Gzj9Mz7+Dq1LvJ3zwbgXSSAGBhGVpMbmcvfTn
NJujNktz5hVtgrNKp+JTZ5MF6hFIsxm53e7wty8/T3ybwXTKIZMpaTPvV14xKHP8OfNv76q0PtK2
T7vhMYqXqJ4egYT2s/nhxNfvjPTgNGq+W6HT+NNVaVHqppXjqSxJpJEOAmlCp0pSp9kwnvi6Tn8j
06h4t0Jb+fNvfIVWwRW8aZZt0QqBhI7EL8yaqXkl2vzw3gnXbwWvrJr0yKHikmSeyogRSHPKLElj
6kW/01pdx38qh8oeYPH4C+bfVpsdal5KtaMT0kgNgXRdwxa7wse++WTUb/wiS3PN9ahKOdE4ZmmO
7XYTIJCmZbvdbk/kvEVC16W5w19X84uK59+R+8LXnBOMCuW34mURSBc1coKIo7FVVWo4/sE5tP5S
wwm6bVXaR6PVboX8kCaNNBFIM7tVkmxnw/qD8Vbjb7U0ZzW7VS5SFVSl9CNl1xwqEUiXY5JG+2g0
3IBnUok2AxBZv2pVlTrdCaJA5jkthaFij0CanNqZpFhZVaqZzRV2K7RNoyZn8ovPKnW9E0Qnsh8H
LATS1Rgem9+KxlMH42Xj71qJ5pjgzlal/R/qbHJ7//Tuxz8/Hn993P9PX798m+DFmhiBNL81CUTm
i71OW7zMl+b2erwEDZPg7tWGFr3vEh36/tffj08HgTTyrhkoQCBhnPT64d2qlD/tKizN7ckeEMQO
q1LmbgWdknRoji47NwLpEpTrUSxx6c+c8QtWojGaJ8FaldbdCmsmuXC4E5000kcgXcLXL99EFity
NlmssfTy+pz5Y60qkflFbDtZ4yf/JQjESxLEEUgQFa/g3ZrjvFSiMVdQbfJbbm1hcF2SqEdeEEjz
Cx9Lnf3f+SN5eX1+//Ru/5e95JAjibNE646YzB+lVpJE3vbIQSBNTmpqKPPy+hxXpcVbDg17CQqS
4O5uhXWf9K194crWM2HwgkC6EI8lKcyA4UyGyOBXOjf4KJC5e3tz6vFwX/ieTkkSOXWKTATSzEQm
hUrro5jjdn+2Tl1u7jBxHVUl8a/fYY9Amtbh59BLSbp1lkj/Hhax8VNhopqczfL0RQ3uViWdkgRH
CCQIiXMosfrvoiqJTMfFV+C+u9glXpXWoxai0RECaU6JT6BUSVrHuTncDn+ezqTFtCqlf7XhDBie
t/CF1qU0s/Of2ERVMkwCkXc4ziKQJuToePDTh89hw0LZgF1UpcHqn5Czs7l4VVooSX4QSFekUJLC
YfXhOedTc4d5Vdozmfs2S3NWp6+Wn48tTJJA6s2AUwik2YgfCe53K2zOVZSNX6cqDX7+m9+ktWY2
V6hKt8ZPSXKBQJpK/kducEkacG2FkVVJ4Rg8HcBl82+Tx7WpSnw7FfkIJHSXnhnjaKw/hrWtSgOO
wZtXoh6sqlI6UClJ+gikeZz9sHUtSfFBcf5NjJpMFol7WHTVe7I7+4jOzr/N3wxxVRpQkhRqKyoR
SJPQOfQruNxciMazdzq4+zOXPpPU4ImvcgN3pk4PKt7Tr4CSJI5AurSGJan+LFGPaWLYCl7zaS6+
I1HNU6ow/4ZhdB0J9WgOBNIMDCedJrsVwlVqOs0pAzY7tH3+Nwk602JXp4u1nxq/SEjjEIHkXuWn
qzgJWk0uY2YHnX3ht3TarXB3/h3ZLRT2hUMZgYRz+m3g7r0TvUlVOrxbYOVTIZ6UrazRmHkPi0wF
LyglSRaB5FuTz1VmEvRYbxk/L7StSjXjH7aBO30JcJMgbFWVOHU0GQLJsTGzeb9KtB//mK/r2l5t
SOeLROOfgU00tq1KNSOBCAIJy3KUBAOurWCrviqV3S/cqpFozr81VYl6NB8Cyat+88uYg9Zb4x95
TaOzVSn+m/nPv04lihku1u2fuoKqVD/+90/vfvzz4/HXx5ofgrYIJJd6pNHL63P4duqYZUCdo/V+
G/CkdivESSDYLU5VpfTdbPN9/+vvxycCSQiB5E/b2XyzNKeQE+PvjnG2Kt3dSB3/WGwk1g8zq9Ld
u9nmj4SSJIVAuq79J//w7kQ9fq9C7O1lVqXEjrX15/QYXr33/95MVnaES0ZVanuwQkmSQiA5Uz+b
p3cr9J6tMsff4+p2ORIXZk3Mg1JLc2kiaXR3k8WtqiS42IiGCCRPKtPIaoutO4kVvPgl0K9Eh2Qb
6saYyzrI7j+8JgJpfqc2cPe+J0X+J9/8Puv7Fbx4U8DiLYeWfyNW5NrbmUkQVyXq0fQIJDfOHsep
fZHI43HopiqNuRNEJ+ujcNcJ1oOATku47p6QiRFIPpz6wFQuzZlXk1XYiW4+Ev0Lszp1NglIjukR
SPNQq0SxCeYR12mkkOs1NvVuaf0mJ+pE/GI9ANx396Oy3v0s/F/9b1w7QRM1H/W2I7mmfRq9H3JP
8Rw5I9mMP75eOCZDQ1KXmM2VK9FMXM99rS5qoKZ5VaIkKaAh+RP6UNtKtNewmlSO0HZjmPdJat2I
saFTMtIjSSw2UpXmQ0OStpkN3X2RqNVkIfJ1TnemOXWU0LAqUZLM0ZB0rZ+NTSUaNgCd8zdWJcn1
9HR3sU6nXlSOhKo0DQJJ2oCluX7azua3lp76cZ1Gi8Uz1tbZehcyqTKWCDZbBJI08xwqLkneZ3Pv
48+czXXm381IyhYbqUreEUi6Pn34LLJipkBn/VCf91NHlSqrEpFmiECSFmZh24m4IAm81wvv4z9F
Z/5dR9LkbrA6jwv5CCRd6yV8fJWDfrP5mOfBexp5r0cNx19clQgzKwSSD7ZVyVciXlnxqZdZ59+1
Ks36ACdDIEmLk8BFVepdL3o/A97rkXefPnzucUnvghW8iUNaGYHkjFVVykkC77O59/HXLHYpzL+9
FxupSvoIJH80q9Kw2VzwsSvwfupo1TUaT1UlhZC+GgJJ3a35d3xVmjsJvNejerbz78hApSrJIpAc
06lKg2fz5o/aexp5r0f7G0z0TovMqkRJGoxAciA9/46sSocj8T6be9d8n3STH+UCVUkNgTQDnao0
TMPHS6DaOgzUYdF4typdLaRtEUg+5My/Y6rSZiTeZ3Pv42++WDd4/hVZbKQqiSCQpjK4KtnO5lcr
hXsis3kng6MxUZXeP7378c+PYSO5MgLJjfz5t3dVsr2Fa0Pe61Enw5JAMFBvVaXvf/1tMp6rIZDm
1LsqffrwWWE2r3mMCuOvITibn3LqbrAjcWFWQwSSJ2fn335Vqcf1XZBvzEUN+v38u3ezNbepSkTU
GATS5Naq1DyWRG5IWlaSvNcj7/LfPIZJQFUaj0Bypmz+bbuC53029z7+MYt1/eZiX4uNa1UinAYg
kC6kSVVaZ3OdTW46IxnA12y+VzB+8ySgKg1DIPlTM/9e8Cu0G97r0UjMwrH3T+9ElqknRiBdUXFV
2szmOtmWORLvaXTBehSYR6Ph7TEv5X/WA0CJ9e7mNT9hOTlBeJ/NvbtsGhmKb48Z/7/ohIZ0afVn
lRyVJNeBarVP2ryarMZf0ygkaPi/xWegukND8qq+JK0/Z8n4sLmezRf/4//04bPr2dDL+DeVKP5z
F+P3jkDCskT14vBTl57NW0VjvXBNI9fBc8j26Q3VpOZZ/f23P5t8k7p+JLfcyiEMRiA51jYJ5liX
ODwS955Srl8RcYnjsPjv8BKMwTkk/GR/VilnNpc6k7Q50+A9jRS2Gtecv2k7m7c6kxTe5OtZovTf
JI2GoSH51mO5LK5KHmfzuCR5HH+s1WKXFcHZPKcSwQqBhGPps0q3/j6f8yZCCfj04fPL63O/Eyen
mA9gVfaElJ0l4i09GIHkXr8kCLOhuyPK9QlRmMfPClG0Gfa6TuXr4SjM5jW7FRTGfzUEEo7F16xb
8j6cUiXJVxrFtzkI/7Fmavzn5g8qfwBd3wnhFq6Pvz6mB7C4OpDCQiDNYUASnF3BsyUVjWmH1efW
4J1WpR6+//X349NBILXawO3l/TMZAgkHDg+Ec6qSSBKIXFwgYV+JDv/O/n9dq1L633aSef7G5D3Q
8IBJ4T18TQQStnK+BrsIV6UwfpFo3MjJocx7ey8CK3iHRt6uaT3HuQi/IZGPQJrEyPk3XZU0k8Bc
Zqc59dSZrOCJ7Prrl0O8ew0RSPhJ/lyjWZXi8StEY04lqqFWlQY84fG7rvnarPkb5uIIpHnUz79n
57VbVckqCXTm5aXoNE/8pJ16LIOr0q2S1PVFH7A0Z3VJdawIJPy/4tlcsyoF46OxuBJVjlOtKrWS
zqG264deLkk+MQJpKlbVZF+VTJLAcC7uvTSXaVhV2idB85d78FEOi3UKCCQsS6PZ3LAqJcbfOxqb
BEDDEZpUpYbjP7s016QkkUYiCKTZFMy/DSevuCrpXAK8h4aVaP961b8iA6pS8+12sgu/GIZAQnuD
o+jutNi2JHm5VsKwqlT53NbvVqiMRuqRDgJpQqfm304T1lqVmv/kjWFrU53OEvWeDbtWpfC91LLx
i3yhlTSSQiBdWu/ZfL3dn/lnvrgkdd2tMOyiBkuf1/rrl28Ft2vqsTQn8nVdVCKQ5qTwndBVv80O
Xecgq6W5Tg+qR1UK+6QzByxSiWI6nxEEBNJ1DTuiXD/2bT//BV/jzRnAsA3c42fDthdmzR//mN0K
Z0sSaSSIQJpWev41Wd9Q/grtMrYSGc6GTVbw4vEnLtwQ/iO+sROQQCChuzgaW1Wlsvn0MKRFvtMa
DDtQ6LfZwXBpLr8kUY80EUgzu1WSzE//VlalJuM3zCGR2bC4Kt0av3L9jYk8/9gjkC7HcLEungV6
nFXKH8nL67NVJKvNhmer0n78UrsV2G7nGoE0OantdnsFVal4uolPYzBhxWqq0tLtThCdKH8cQCBd
i+HB461oPFWVCsa/X5rL36ncXPph2h7a51SlMP7DSqRTTRIjIY3EEUjzW5NAZL441GMDXmJ6NbnR
gP5smFOV9B8F/CKQME56/fBuVcoM1JzdCuPj2dHN3/ZVKfMskXhJIkr1EUiXIH4mKXarKuXMdKdO
zg8uSb5u/haewziHvF+43cv7/+IIpEvQOeGcE41rVcq8TprUF4kOZV4hQmf86zFBeG43p44SxEsS
xBFIV+GoJAVxVTqcWSpzaNgT4mixbr80N8eE7uudf2UE0vzUjhPzkyB8W+jWtRWkHlSC/mJd+hTR
uiMm86fpVJMwEv3nHysCaXLx1OCuJC3/1oswY64reA3vb9v7CRF/wu/ubFzHb3JndFwNgQQDZ5Mg
RNGnD58XjcsBZDK/TeItxddWyLysg05Joh75QiDNbD8pOC1J4VH0uNqQxyekxqkvex0+M23vYdGV
/tfvsEEgTUv8c5hIgs3pisNMXeSrklTOFVSi9PjvruDplCQ4QiBBSHzcvU5n+4PxtlWpR0kSSaPe
lz3tdw+LevHZL6LRCwJpTolPoM4iVXxNo/Anh2O+dTDupSrd1WO6rHxm8t8hiapkmAQi73CcRSBN
yNfxYOZoDw/GW1WltiFtOBs2qUQF41euSgslyQ8CaTY5HzzzkrRWolMjcVGVrJ5Y82fg8NUxSQLq
kV8EEsbJv7bC3bPlS+uqZBXSlfN187NElU+CeVW6NX5KkgsE0lTyP3KD59/EJFUwEtmqNOwp7bRb
ocn4N6/OGlEDOLpKEw4RSPMQPAAMM1GngTWvSpUhPSaNzJfmMplUpfTXYClJ+gik6+pakuLJKOec
VtlMIVuVmuu9gXvpEKjxqzOgJHHqaAIE0iREDv2Kr8AdH9uefSwNq1JxSJf9q7uPdEAOrb+o93eV
usocPyVJHIE0g+LPWKuSVH9HospLvNhWpR6z+TQNb31vkAS4i0BClYbnCeqvg5muSkve/G67J35Y
Jdr80h6/bn/djU5nlU6Nn5KkjEByr/LTVTb/yt6kdczVhmKtfqZVJWr+nCTeGz3uYcGpo5kQSL6N
P9brt3UqRGO4KV/lj6pcwcsP6ZrZMAzPpBLFY2j1exM5dPht2cO/OQYlSRaBhKz5V7YSJXS92lAl
2xxaNbldUEG6tLqHhfnriLYIJMcGHOUNzqHm52+Kq1LOSIr34y3VmziaqHyq698blSt4NeOnJGki
kLxq+3Haz78miyqd9gf3qEpn/6FIJVoVP/DmxyjmK3jQQSDhJ+ZLc502uRVUpVYjOfz5Tg/PC4ad
808KqlL9S0NJEkQgudTjg/Ty+jz9inyrqpTzl9UqUezUg41ra9fpO78qTf9GvSwCyZ+2aWReifa6
fhPoVFU6PNOTHptyDgX5z+34lbQe+8LTv46SJIVAuqjDHLL9TujIqSG/Km32oSX2Sed/kUh8EjQ/
RklXJerRxAgkZ+rnMhcnkAdEY2ZV+hjdYX052ietX4liiWe1Rw4Vv11vvTo9rgArfnxwKQSSJzWf
nMzpxrYkjZdTlT59+Bz+fPPMuLvc3K1XVvYYhQ14V0MgTc58+aXYsGi8W5WCdbHOVyVKcPHeiF+d
Tu8HSpIOAsmNs5+Z4kNLk5JkPiPcvTDr2pMqnxmrRxq/pmNqR8NHGtLo5fW5yU+DLALJh/zPtovD
3kyDozFdlVz3oTB+7++NflXm/dO7H//8ePz1sflPxikE0jwaHvZOMAUXm++8RdcbyY+xvhv7vTrf
//r78YlAMkYgOZCeSrwf9qaZRGOrS3+aCw+hyRVUDW3eANO8OtgjkNTdSqOZckjz4L3HlzTHPNL4
vWHVdHs/0uavDlsbFBBI/ow5NtRZtbMdybpG5MJ8vSH90s+3vnpxBJK0+JBtpkrky+Dr2RS49d4Q
OaQoljP+hq8OJckcgaRr/WxYHQNSkmKCB+PpYxTzZ2wkwVcHBQgkaVf4jDk6Jq08GG/4SMXfGE0e
6dlAbVKVKEm2frEeAG769OHz+6d3tp+NzdUKDOmMJMxZJieWwu8NM2b6jeG9HhWP3/DVQT0akqi3
t7dlWR4eHlxPK7Maf1bpVCXynkaJS6rnqHx1KEmGaEjS3t7efv/tT9tm0LWanPrk65SkYMDBeH4l
msl6laYaVCWPaEjqqErK8g/GT0VvzY5K83pUWS8ajr+4KlGSrBBIPry9vT08PCxGV5hW2OS2jkRw
pmi4xavy54i8TMV6jJ8NeI4QSG5QlQLNC+HUb8CLfw7aKnh1KEkmCCRnrKpStq7gEwAAB3lJREFU
j5I03wf+1sF44pE2PHinHqVRlfQRSP6sVWnxf4+4Ajrrh4cyv8vcvBKJPCfFBxm2N2O89ZfnO2YS
RyB5NX4FTzwJpNya+FiaU0BVkkUg+Wa72cGKl2iML8zadQZ08WwkyN5hhJI0GN9Dcu/t7S18XWnA
72r4TaCrfc77fZFIJ43KXlPbS7n7upr79AikSSh8hXYktS/J7q1faF16fkmz8qIGWO69OiTWSATS
PMZUJf0ksHXr2gqdDsabXNTAkEi9oyqJ4BzSbK5zVknqTFLmboW2p9N1Hn4ZtfHfenU4kzQMgTSh
3hvw6pNgpo93OmD2j7TVhVnVZvM5XlP9mzHOjUCa1hWqkmFJqtzAzc5jtUCN7V8dStIYBNLM+lUl
qeWykRp+kajmYNz7k68/fqqSCQJpfnNXpWHR2KnQFFQl/dl8GvGrw66HAQikS+hRlYqTwNdRZ00l
yr9ETf5f1nRq8L4CdYJXxxEC6UJmrUqdStLgczyZVcnXbL7ndPzUozEIpGtpe2HWKc8kGV5u7u7B
+HzPtiPvn97xDbzeCKQrMry1UqeljyY70cN/mK/MTLwBz2mgrjkUPjjoh0C6riYreJq3cD2l0+xf
87QcXvpTdjbPfKSy408IUUQODUMgXVqTqiRyC9ezJUmnEt0Sr+B5nM39ohJZIZAw7WaHW3wtiM1x
Ot1FoJJD5ggkLEtdVcqvJr0X99IjGVmJ2j7S0EH9Lo3qpxFLcyIIJPxnyqqkvzSXts7mmpsd/Mbk
QiXSQyDhJ2VVSWf/dzwSwem7kscvaYq8MTaoRJoIJBzwXpWmiaLD2VyzKh1SSyMqkTgCCcfOVqW7
JWnAcX2Yps3rWqtHmngUIlXJfACZyCEvCCSkuKhK3s8SFROvSgr1iKU5Xwgk3JFflQZXk0QOud6T
FmQ+kyJVac82jahEThFIyCJVlXJqgcjXdcucnc3VqtLXL9+snnwqkWsEEnLlVKVbJanJIfyppTmr
kmRVVsZXpcTvGn80QCWaA4GEc0yqUtnhv9OSVLPYpVCVRi7WkUOTIZBwWvoeFg3PJF1wt0L9U2d7
VmlYGrE0NyUCCYX63cOiVQ6Z7/82pFCVeqASzY1AQpXDFbxNEuQfrXufQyt7Sdv4PLyHRSuHj7Rr
/FOJroBAQq36qtRvac5RSeo0zmEreJ3GTyW6FAIJbWyqUmYSeK9ErfTeJ+1uBY8cuiYCCc3cqkr7
SXDkbgUXJWnAhsCuVanhM8zS3JURSGgsrkrrcW4w96654rl+ZF42qUqbR9pk/FQiLAQSeoj3hQe2
60XKJWn8wKSuNkQOIUYgoZe1KonMfYi1OqtUHKgszWGPQEJH/b6rdFbvklQWura9rb4qFYyfSoQE
AgndSV2YVYfIKuLZqlRzqmwhh5BEIGEEhaqkfCbJVllVynwyqUTIRyBhHKrSSjAaT1Wlu+Mnh1CA
QMJQtlWpU0nq1C3Gu1uVch4pS3MoRiDBAFVJ2d2qdBioVCLUI5BgI30Pi37MzyTJ1qNYoirtx08l
QisEEiwpbHYYyUUardJViUqE5ggk2Bu8gte2JM39td/4HhbLsnz68HmhEqEbAgkSrlCVfNWjWJy4
5BD6IZAgZFhVGn8myW8asTSHYQgkaLlCVfKCpTkMRiBB0YCq1KQkZZ5A8lWPqESwQiBB1DRVyUsa
kUMwRyBBWteqZP6dJBEszUEEgQR1rquScuBRiaCGQIIPnarSx+9/FH+R6O4/lE0jKhE0EUhwo1NV
+vThs2ZsNEclgjgCCc54uTCrTj0ih+AFgQR/2lalsq0N6fU6kTRiaQ6+EEjwyktVGo9KBKcIJDjW
6h4Wbfd/G9YjKhFcI5DgntS+cJM0ohJhDgQSJlG5gneqJInccoIcwmQIJMzDvCoNq0cszWFKBBJm
U1yVKs8kDUgjKhHmRiBhQuOr0tcv3/r9LnIIF0EgYVoFVSmnJB2eQOp0uQeW5nApBBJmNqYqNV+s
oxLhmggkzO9UVTp7JqltGlGJcGUEEi7BfANeGpUIWAgkXEpmVUqUpM0JpMp6RA4BMQIJ19KwKlXu
EV/IIeBnBBKu6G5V6nR3cyoRkEAg4aIqL8x6Nq6oRMBdBBIuLbGCty9J6wmk/DSiEgH5CCSg/a2V
yCGgAIEELMuNqnR4Jildj1iaA4oRSMB/ElUprNfdSiMqEVCPQAJ+sqlKoSQl/j6VCGiFQAIOHFal
uB5RiYDmHvg4AQmhKoX4idsSHxygOQIJuCNUpYDPC9APgQQAkPCL9QAAAFgWAgkAIIJAAgBIIJAA
ABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA
ABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA
ABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA
ABIIJACABAIJACCBQAIASCCQAAASCCQAgIT/A/atmRzsFwPsAAAAAElFTkSuQmCC
"
>
</div>

</div>

</div>
</div>

</div>
<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>For 3-D linear element, <code>mg</code> wins at an even smaller size $3.6\times 10^4$. Again <code>amg</code> is 3-4 times slower than <code>mg</code>.</p>

</div>
</div>
</div>
    </div>
  </div>
</body>

 


</html>
